Cross species correlation - pairwise genes (updated orthology)

Notebook for cell type identity correlations using pairwise genes. Subsetting will be needed since for some of the species we only used the first cell type name (which is what we’re more interested in).

General Setup

Setup chunk

Setup reticulate

knitr::opts_chunk$set(fig.width = 8)
knitr::opts_knit$set(root.dir = normalizePath(".."))
knitr::opts_knit$get("root.dir")

[1] "/nas/groups/treutlein/USERS/tomasgomes/projects/pallium_evo"

Load libraries

library(reticulate)
knitr::knit_engines$set(python = reticulate::eng_python)
py_available(initialize = FALSE)

[1] FALSE

use_python(Sys.which("python"))
py_config()

python:         /home/tpires/bin/miniconda3/bin/python
libpython:      /home/tpires/bin/miniconda3/lib/libpython3.8.so
pythonhome:     /home/tpires/bin/miniconda3:/home/tpires/bin/miniconda3
version:        3.8.3 (default, May 19 2020, 18:47:26)  [GCC 7.3.0]
numpy:          /home/tpires/bin/miniconda3/lib/python3.8/site-packages/numpy
numpy_version:  1.18.5

NOTE: Python version was forced by RETICULATE_PYTHON

Data

Load datasets

library(Seurat)

Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     
Registered S3 method overwritten by 'htmlwidgets':
  method           from         
  print.htmlwidget tools:rstudio
Attaching SeuratObject

library(simspec)
library(harmony)

Loading required package: Rcpp

library(dplyr)

Attaching package: ‘dplyr’

The following objects are masked from ‘package:stats’:

    filter, lag

The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union

library(tidyverse)

Registered S3 methods overwritten by 'dbplyr':
  method         from
  print.tbl_lazy     
  print.tbl_sql      
Registered S3 method overwritten by 'cli':
  method     from         
  print.boxx spatstat.geom
── Attaching packages ─────────────────────────────────────────────────────────── tidyverse 1.3.0 ──
✓ ggplot2 3.3.3     ✓ purrr   0.3.4
✓ tibble  3.1.2     ✓ stringr 1.4.0
✓ tidyr   1.1.3     ✓ forcats 0.5.0
✓ readr   1.4.0     
── Conflicts ────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
x dplyr::filter() masks stats::filter()
x dplyr::lag()    masks stats::lag()

library(ggplot2)
library(uwot)

Loading required package: Matrix

Attaching package: ‘Matrix’

The following objects are masked from ‘package:tidyr’:

    expand, pack, unpack

Load metadatas

lizard_all = readRDS("data/expression/lizard_all_v3.RDS")
turtle_all = readRDS("data/expression/turtle_all_v3.RDS")
axolotl_all = readRDS(file = "data/expression/axolotl/palliumres07_harmony.RDS")

Add metadata

lizard_meta = read.csv("data/annotations/lizard_all_umeta.csv", header = T, row.names = 1)
turtle_meta = read.csv("data/annotations/turtle_all_umeta.csv", header = T, row.names = 1)
axolotl_meta = read.csv("data/annotations/axolotl_all_umeta.csv", header = T, row.names = 1)

Filter data (doublets and odd cells)

lizard_all = AddMetaData(lizard_all, metadata = lizard_meta)
turtle_all = AddMetaData(turtle_all, metadata = turtle_meta)
axolotl_all = AddMetaData(axolotl_all, metadata = axolotl_meta)

Put Seurats in list

turtle_all = turtle_all[,turtle_all$classes!="doublets"]
axolotl_all = axolotl_all[,!axolotl_all$cellclusters %in% c("GABAergic neuron", "ETV1+ neuron", "neuron")]

Ortholog matching

Load necessary orthology data

srat_list = list("paintedturtle" = turtle_all, "axolotl" = axolotl_all, "lizard" = lizard_all)

Function for ortholog matching

ort_alt_ens_list = readRDS(file = "data/eggNOG/ort_alt_ens_list.RDS")
ort_alt_egg_list = readRDS(file = "data/eggNOG/ort_alt_egg_list.RDS")

ort_list_use = c(ort_alt_ens_list, 
                 ort_alt_egg_list[which(grepl("axolotl", names(ort_alt_egg_list)) |
                                          grepl("lizard", names(ort_alt_egg_list)))])
bfmats = ort_alt_egg_list[which(grepl("zebrafinch", names(ort_alt_egg_list)))]
names(bfmats) = gsub("zebrafinch", "bengalesefinch", names(bfmats))
ort_list_use = c(ort_list_use, bfmats)

Function for subsetting Seurats

# get orthologs for all species pairs
ortMatch = function(ort_l, sp_vec, ort_scope = "ortholog_one2one"){
  sp_p_mat = combn(sp_vec, 2)
  for(i in 1:ncol(sp_p_mat)){
    sp_pair = sp_p_mat[,i]
    ort_tab = ort_l[[paste0(sp_pair, collapse = ".vs.")]]
    if(!is.null(ort_scope)){
      ort_tab = ort_tab[ort_tab$homology.type==ort_scope,]
    }

    colnames(ort_tab) = gsub("genename", "", colnames(ort_tab))
    if(i==1){
      ort_all = unique(ort_tab[,c(2,4)])
    } else{
      ort_all = unique(merge(ort_all, unique(ort_tab[,c(2,4)]), 
                             by = intersect(colnames(ort_all), colnames(ort_tab))))
    }
  }
  
  return(ort_all)
}

Get ortholog-matched Seurats

seuratOrthologs = function(s_l, o_all){
  # only orthologs present in the Seurat objects
  for(sp in names(s_l)){
    o_all = o_all[o_all[,sp] %in% rownames(s_l[[sp]]),]
  }
  
  # get a unique name for the genes
  allg = unlist(o_all)
  jointnames = c()
  rep_n = list()
  for(i in 1:nrow(o_all)){
    nname = unlist(lapply(o_all[i,], function(x) sum(allg==x)))
    g = unlist(o_all[i,names(nname[nname==1])])
    # if there are acceptable unique names
    if(any(nname==1) & (!is.null(g) & !all(grepl("..", g, fixed = T)))){ # no weird axololt annot
      g = g[!grepl("..", g, fixed = T)] # no weird axololt annot
      g = c(g[!grepl("LOC", g)], g[grepl("LOC", g)])[1] #prioritize non-LOC
      jointnames = c(jointnames, g)
      rep_n[[g]] = 1
    } else{ # we'll add a number to those names
      g = unlist(o_all[i,])
      g = g[!grepl("..", g, fixed = T)] # no weird axololt annot
      g = c(g[!grepl("LOC", g)], g[grepl("LOC", g)])[1] #prioritize non-LOC
      if(g %in% names(rep_n)){
        gn = rep_n[[g]]
        rep_n[[g]] = gn+1
        jointnames = c(jointnames, paste0(g, "-", gn+1))
      } else{
        rep_n[[g]] = 1
        jointnames = c(jointnames, paste0(g, "-", 1))
      }
    }
  }
  
  # filter the Seurats by the orthologs
  ## this requires redoing the Seurat with new names, since some genes will be duplicated
  for(sp in names(s_l)){
    cnts = s_l[[sp]]@assays$RNA@counts[o_all[,sp],]
    rownames(cnts) = jointnames
    m = s_l[[sp]]@meta.data
    s_l[[sp]] = CreateSeuratObject(cnts, project = sp, meta.data = m)
  }
  
  return(s_l)
}

Integration

Normalise data

for(sp in names(srat_ort)){
  print(sp)
  srat_ort[[sp]] = suppressWarnings(SCTransform(srat_ort[[sp]], do.correct.umi = T, verbose = F,
                                                seed.use = 1, vars.to.regress = "nCount_RNA",
                                                variable.features.rv.th = 1, return.only.var.genes = F,
                                                variable.features.n = NULL))
  srat_ort1[[sp]] = suppressWarnings(SCTransform(srat_ort1[[sp]], do.correct.umi = T, verbose = F,
                                                 seed.use = 1, vars.to.regress = "nCount_RNA",
                                                 variable.features.rv.th = 1, return.only.var.genes = F,
                                                 variable.features.n = NULL))
}

Seurats with all merged

for(sp in names(srat_ort)){
  print(sp)
  srat_ort[[sp]] = suppressWarnings(SCTransform(srat_ort[[sp]], do.correct.umi = T, verbose = F,
                                                seed.use = 1, vars.to.regress = "nCount_RNA",
                                                variable.features.rv.th = 1, return.only.var.genes = F,
                                                variable.features.n = NULL))
  srat_ort1[[sp]] = suppressWarnings(SCTransform(srat_ort1[[sp]], do.correct.umi = T, verbose = F,
                                                 seed.use = 1, vars.to.regress = "nCount_RNA",
                                                 variable.features.rv.th = 1, return.only.var.genes = F,
                                                 variable.features.n = NULL))
}

[1] "paintedturtle"
[1] "axolotl"
[1] "lizard"

Seurat rPCA

Integration with all ortholog

srat_ort_all = Reduce(merge, srat_ort)
srat_ort_all = AddMetaData(srat_ort_all, 
                           metadata = unlist(lapply(names(srat_ort), 
                                                    function(x) rep(x, ncol(srat_ort[[x]])))),
                           col.name = "species")

srat_ort1_all = Reduce(merge, srat_ort1)
srat_ort1_all = AddMetaData(srat_ort1_all, 
                           metadata = unlist(lapply(names(srat_ort1), 
                                                    function(x) rep(x, ncol(srat_ort1[[x]])))),
                           col.name = "species")

Integration with one2one ortholog

feat_list = list("all" = rownames(srat_ort1_all@assays$SCT@scale.data), 
                 "allvar" = unique(unlist(lapply(srat_ort1, VariableFeatures))), 
                 "commonvar" = Reduce(intersect, lapply(srat_ort1, VariableFeatures)))
varindat = lapply(srat_ort1, function(x) feat_list$allvar %in% rownames(x))
feat_list$allvar = feat_list$allvar[rowSums(Reduce(cbind, varindat))==length(srat_ort1)]
srat_rpca_one_l = list()
for(n in names(feat_list)){
  print(n)
  integr_feat = feat_list[[n]]
  
  # calculate possible missing Pearson residuals for SCTransform
  srat_ort1_int = PrepSCTIntegration(srat_ort1, anchor.features = integr_feat, verbose = T)
  
  srat_ort1_int <- lapply(X = srat_ort1_int, FUN = function(x) {
    x <- ScaleData(x, features = integr_feat, verbose = FALSE)
    x <- RunPCA(x, features = integr_feat, verbose = FALSE)
  })
  
  # finding the anchors for integration
  all_cell_anchors = FindIntegrationAnchors(srat_ort1_int, normalization.method = "SCT", 
                                            anchor.features = integr_feat, reduction = "rpca",
                                            assay = rep("SCT", length(srat_ort1_int)), 
                                            dims = 1:50, verbose = T)
  
  # actual data integration
  allgenes = rownames(srat_ort1_int[[1]]@assays$RNA@counts)
  all_cell_integr = IntegrateData(all_cell_anchors, normalization.method = "SCT", dims = 1:50,
                                  verbose = T, features.to.integrate = allgenes)
  
  # run PCA and UMAP to see how it looks
  all_cell_integr = RunPCA(all_cell_integr, verbose = F)
  all_cell_integr = RunUMAP(all_cell_integr, dims = 1:30)
  
  all_cell_integr = AddMetaData(all_cell_integr, 
                             metadata = unlist(lapply(names(srat_ort1), 
                                                      function(x) rep(x, ncol(srat_ort1[[x]])))),
                             col.name = "species")
  srat_rpca_one_l[[n]] = all_cell_integr
}

Harmony

Integration with all ortholog

feat_list = list("all" = rownames(srat_ort_all@assays$SCT@scale.data), 
                 "allvar" = unique(unlist(lapply(srat_ort, VariableFeatures))), 
                 "commonvar" = Reduce(intersect, lapply(srat_ort, VariableFeatures)))
varindat = lapply(srat_ort, function(x) feat_list$allvar %in% rownames(x))
feat_list$allvar = feat_list$allvar[rowSums(Reduce(cbind, varindat))==length(srat_ort)]
srat_harm_all_l = list()
for(n in names(feat_list)){
  print(n)
  # Renormalise
  integr_feat = feat_list[[n]]
  VariableFeatures(srat_ort_all) = integr_feat
  srat_ort_all_int = ScaleData(srat_ort_all, features = integr_feat, use.umi = T,
                               do.scale = F, verbose = FALSE)
  srat_ort_all_int = RunPCA(srat_ort_all_int, verbose = FALSE, assay = "SCT", npcs = 50,
                            features = integr_feat)
  
  # Run Harmony
  srat_ort_all_int = RunHarmony(srat_ort_all_int, "species", tau = 30, 
                                plot_convergence = F, assay.use = "SCT")
  
  # Run UMAP on Harmony dimensions
  srat_ort_all_int = RunUMAP(srat_ort_all_int, reduction = "harmony", dims = 1:10)
  
  srat_harm_all_l[[n]] = srat_ort_all_int
}

Integration with one2one ortholog

feat_list = list("all" = rownames(srat_ort1_all@assays$SCT@scale.data), 
                 "allvar" = unique(unlist(lapply(srat_ort1, VariableFeatures))), 
                 "commonvar" = Reduce(intersect, lapply(srat_ort1, VariableFeatures)))
varindat = lapply(srat_ort1, function(x) feat_list$allvar %in% rownames(x))
feat_list$allvar = feat_list$allvar[rowSums(Reduce(cbind, varindat))==length(srat_ort1)]
srat_harm_one_l = list()
for(n in names(feat_list)){
  print(n)
  # Renormalise
  integr_feat = feat_list[[n]]
  VariableFeatures(srat_ort1_all) = integr_feat
  srat_ort1_all_int = ScaleData(srat_ort1_all, features = integr_feat, use.umi = T,
                               do.scale = F, verbose = FALSE)
  srat_ort1_all_int = RunPCA(srat_ort1_all_int, verbose = FALSE, assay = "SCT", npcs = 50,
                            features = integr_feat)
  
  # Run Harmony
  srat_ort1_all_int = RunHarmony(srat_ort1_all_int, "species", tau = 30, 
                                 plot_convergence = F, assay.use = "SCT")
  
  # Run UMAP on Harmony dimensions
  srat_ort1_all_int = RunUMAP(srat_ort1_all_int, reduction = "harmony", dims = 1:10)
  
  srat_harm_one_l[[n]] = srat_ort1_all_int
}

CSS

Integration with all ortholog

feat_list = list("all" = rownames(srat_ort_all@assays$SCT@scale.data), 
                 "allvar" = unique(unlist(lapply(srat_ort, VariableFeatures))), 
                 "commonvar" = Reduce(intersect, lapply(srat_ort, VariableFeatures)))
varindat = lapply(srat_ort, function(x) feat_list$allvar %in% rownames(x))
feat_list$allvar = feat_list$allvar[rowSums(Reduce(cbind, varindat))==length(srat_ort)]
srat_css_all_l = list()
for(n in names(feat_list)){
  print(n)
  integr_feat = feat_list[[n]]
  VariableFeatures(srat_ort_all) = integr_feat
  srat_ort_all_int = ScaleData(srat_ort_all, features = integr_feat, use.umi = T,
                               do.scale = F, verbose = FALSE)
  srat_ort_all_int = RunPCA(srat_ort_all_int, verbose = FALSE, assay = "SCT", npcs = 50,
                            features = integr_feat)
  
  # Run CSS
  srat_ort_all_int = cluster_sim_spectrum(srat_ort_all_int, label_tag="species", redo_pca = T, 
                                          dims_use = 1:30, var_genes = integr_feat)
  
  # Run UMAP on CSS dimensions
  srat_ort_all_int = RunUMAP(srat_ort_all_int, reduction="css", reduction.name="umap_css",
                             dims = 1:ncol(Embeddings(srat_ort_all_int,"css")), reduction.key="UMAPCSS_")
  srat_ort_all_int = SetIdent(srat_ort_all_int, value = "species")
  
  srat_css_all_l[[n]] = srat_ort_all_int
}

Integration with one2one ortholog

feat_list = list("all" = rownames(srat_ort_all@assays$SCT@scale.data), 
                 "allvar" = unique(unlist(lapply(srat_ort, VariableFeatures))), 
                 "commonvar" = Reduce(intersect, lapply(srat_ort, VariableFeatures)))
varindat = lapply(srat_ort, function(x) feat_list$allvar %in% rownames(x))
feat_list$allvar = feat_list$allvar[rowSums(Reduce(cbind, varindat))==length(srat_ort)]
srat_css_all_l = list()
for(n in names(feat_list)){
  print(n)
  integr_feat = feat_list[[n]]
  VariableFeatures(srat_ort_all) = integr_feat
  srat_ort_all_int = ScaleData(srat_ort_all, features = integr_feat, use.umi = T,
                               do.scale = F, verbose = FALSE)
  srat_ort_all_int = RunPCA(srat_ort_all_int, verbose = FALSE, assay = "SCT", npcs = 50,
                            features = integr_feat)
  
  # Run CSS
  srat_ort_all_int = cluster_sim_spectrum(srat_ort_all_int, label_tag="species", redo_pca = T, 
                                          dims_use = 1:30, var_genes = integr_feat)
  
  # Run UMAP on CSS dimensions
  srat_ort_all_int = RunUMAP(srat_ort_all_int, reduction="css", reduction.name="umap_css",
                             dims = 1:ncol(Embeddings(srat_ort_all_int,"css")), reduction.key="UMAPCSS_")
  srat_ort_all_int = SetIdent(srat_ort_all_int, value = "species")
  
  srat_css_all_l[[n]] = srat_ort_all_int
}

[1] "all"
Start to do clustering for each sample...
>> Redoing truncated PCA on sample axolotl...
found nearest neighbors.
revoke Seurat to compute SNN.
done.
done. returning result...
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 16415
Number of edges: 611282

Running Louvain algorithm...

0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

Maximum modularity in 10 random starts: 0.9237
Number of communities: 20
Elapsed time: 2 seconds
>> Done clustering of sample axolotl.
>> Redoing truncated PCA on sample lizard...
found nearest neighbors.
revoke Seurat to compute SNN.
done.
done. returning result...
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 4187
Number of edges: 160419

Running Louvain algorithm...

0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

Maximum modularity in 10 random starts: 0.8918
Number of communities: 11
Elapsed time: 0 seconds
>> Done clustering of sample lizard.
>> Redoing truncated PCA on sample paintedturtle...
found nearest neighbors.
revoke Seurat to compute SNN.
done.
done. returning result...
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 18664
Number of edges: 619386

Running Louvain algorithm...

0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

Maximum modularity in 10 random starts: 0.8942
Number of communities: 15
Elapsed time: 2 seconds

1 singletons identified. 14 final clusters.

>> Done clustering of sample paintedturtle.
Finished clustering.
Obtained average profiles of clusters.
Start to calculate standardized similarities to clusters...
Doing z-transformation...
Done.

23:28:48 UMAP embedding parameters a = 0.9922 b = 1.112
23:28:48 Read 39266 rows and found 45 numeric columns
23:28:48 Using Annoy for neighbor search, n_neighbors = 30
23:28:48 Building Annoy index with metric = cosine, n_trees = 50
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
23:28:57 Writing NN index file to temp file /tmp/RtmpyIf2G8/file2bdbd14969e74
23:28:57 Searching Annoy index using 1 thread, search_k = 3000
23:29:10 Annoy recall = 100%
23:29:13 Commencing smooth kNN distance calibration using 1 thread
23:29:17 Initializing from normalized Laplacian + noise
23:29:21 Commencing optimization for 200 epochs, with 1762176 positive edges
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
23:30:11 Optimization finished

[1] "allvar"
Start to do clustering for each sample...
>> Redoing truncated PCA on sample axolotl...
found nearest neighbors.
revoke Seurat to compute SNN.
done.
done. returning result...
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 16415
Number of edges: 613106

Running Louvain algorithm...

0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

Maximum modularity in 10 random starts: 0.9233
Number of communities: 19
Elapsed time: 2 seconds
>> Done clustering of sample axolotl.
>> Redoing truncated PCA on sample lizard...
found nearest neighbors.
revoke Seurat to compute SNN.
done.
done. returning result...
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 4187
Number of edges: 160511

Running Louvain algorithm...

0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

Maximum modularity in 10 random starts: 0.8930
Number of communities: 11
Elapsed time: 0 seconds
>> Done clustering of sample lizard.
>> Redoing truncated PCA on sample paintedturtle...
found nearest neighbors.
revoke Seurat to compute SNN.
done.
done. returning result...
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 18664
Number of edges: 615768

Running Louvain algorithm...

0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

Maximum modularity in 10 random starts: 0.8890
Number of communities: 14
Elapsed time: 2 seconds

1 singletons identified. 13 final clusters.

>> Done clustering of sample paintedturtle.
Finished clustering.
Obtained average profiles of clusters.
Start to calculate standardized similarities to clusters...
Doing z-transformation...
Done.

23:40:23 UMAP embedding parameters a = 0.9922 b = 1.112
23:40:23 Read 39266 rows and found 43 numeric columns
23:40:23 Using Annoy for neighbor search, n_neighbors = 30
23:40:23 Building Annoy index with metric = cosine, n_trees = 50
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
23:40:29 Writing NN index file to temp file /tmp/RtmpyIf2G8/file2bdbd16b007fa
23:40:29 Searching Annoy index using 1 thread, search_k = 3000
23:40:42 Annoy recall = 100%
23:40:44 Commencing smooth kNN distance calibration using 1 thread
23:40:48 Initializing from normalized Laplacian + noise
23:40:52 Commencing optimization for 200 epochs, with 1764938 positive edges
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
23:41:42 Optimization finished

[1] "commonvar"
Start to do clustering for each sample...
>> Redoing truncated PCA on sample axolotl...
found nearest neighbors.
revoke Seurat to compute SNN.
done.
done. returning result...
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 16415
Number of edges: 612295

Running Louvain algorithm...

0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

Maximum modularity in 10 random starts: 0.9173
Number of communities: 16
Elapsed time: 2 seconds
>> Done clustering of sample axolotl.
>> Redoing truncated PCA on sample lizard...
found nearest neighbors.
revoke Seurat to compute SNN.
done.
done. returning result...
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 4187
Number of edges: 169158

Running Louvain algorithm...

0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

Maximum modularity in 10 random starts: 0.8871
Number of communities: 10
Elapsed time: 0 seconds
>> Done clustering of sample lizard.
>> Redoing truncated PCA on sample paintedturtle...
found nearest neighbors.
revoke Seurat to compute SNN.
done.
done. returning result...
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 18664
Number of edges: 642681

Running Louvain algorithm...

0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

Maximum modularity in 10 random starts: 0.8838
Number of communities: 14
Elapsed time: 2 seconds
>> Done clustering of sample paintedturtle.
Finished clustering.
Obtained average profiles of clusters.
Start to calculate standardized similarities to clusters...
Doing z-transformation...
Done.

23:44:05 UMAP embedding parameters a = 0.9922 b = 1.112
23:44:05 Read 39266 rows and found 40 numeric columns
23:44:05 Using Annoy for neighbor search, n_neighbors = 30
23:44:05 Building Annoy index with metric = cosine, n_trees = 50
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
23:44:11 Writing NN index file to temp file /tmp/RtmpyIf2G8/file2bdbd148aac04
23:44:11 Searching Annoy index using 1 thread, search_k = 3000
23:44:24 Annoy recall = 100%
23:44:25 Commencing smooth kNN distance calibration using 1 thread
23:44:29 Initializing from normalized Laplacian + noise
23:44:33 Commencing optimization for 200 epochs, with 1796060 positive edges
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
23:45:23 Optimization finished

Combined

rPCA + Harmony, all ortholog

feat_list = list("all" = rownames(srat_ort1_all@assays$SCT@scale.data), 
                 "allvar" = unique(unlist(lapply(srat_ort1, VariableFeatures))), 
                 "commonvar" = Reduce(intersect, lapply(srat_ort1, VariableFeatures)))
varindat = lapply(srat_ort1, function(x) feat_list$allvar %in% rownames(x))
feat_list$allvar = feat_list$allvar[rowSums(Reduce(cbind, varindat))==length(srat_ort1)]
srat_css_one_l = list()
for(n in names(feat_list)){
  print(n)
  integr_feat = feat_list[[n]]
  VariableFeatures(srat_ort1_all) = integr_feat
  srat_ort1_all_int = ScaleData(srat_ort1_all, features = integr_feat, use.umi = T,
                               do.scale = F, verbose = FALSE)
  srat_ort1_all_int = RunPCA(srat_ort1_all_int, verbose = FALSE, assay = "SCT", npcs = 50,
                            features = integr_feat)
  
  # Run CSS
  srat_ort1_all_int = cluster_sim_spectrum(srat_ort1_all_int, label_tag="species", redo_pca = T, 
                                          dims_use = 1:30, var_genes = integr_feat)
  
  # Run UMAP on CSS dimensions
  srat_ort1_all_int = RunUMAP(srat_ort1_all_int, reduction="css", reduction.name="umap_css",
                             dims = 1:ncol(Embeddings(srat_ort1_all_int,"css")),
                             reduction.key="UMAPCSS_")
  srat_ort1_all_int = SetIdent(srat_ort1_all_int, value = "species")
  
  srat_css_one_l[[n]] = srat_ort1_all_int
}

[1] "all"
Start to do clustering for each sample...
>> Redoing truncated PCA on sample axolotl...
found nearest neighbors.
revoke Seurat to compute SNN.
done.
done. returning result...
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 16415
Number of edges: 640464

Running Louvain algorithm...

0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

Maximum modularity in 10 random starts: 0.9261
Number of communities: 20
Elapsed time: 2 seconds
>> Done clustering of sample axolotl.
>> Redoing truncated PCA on sample lizard...
found nearest neighbors.
revoke Seurat to compute SNN.
done.
done. returning result...
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 4187
Number of edges: 163557

Running Louvain algorithm...

0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

Maximum modularity in 10 random starts: 0.8957
Number of communities: 11
Elapsed time: 0 seconds
>> Done clustering of sample lizard.
>> Redoing truncated PCA on sample paintedturtle...
found nearest neighbors.
revoke Seurat to compute SNN.
done.
done. returning result...
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 18664
Number of edges: 626082

Running Louvain algorithm...

0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

Maximum modularity in 10 random starts: 0.8882
Number of communities: 13
Elapsed time: 3 seconds
>> Done clustering of sample paintedturtle.
Finished clustering.
Obtained average profiles of clusters.
Start to calculate standardized similarities to clusters...
Doing z-transformation...
Done.

23:55:34 UMAP embedding parameters a = 0.9922 b = 1.112
23:55:34 Read 39266 rows and found 44 numeric columns
23:55:34 Using Annoy for neighbor search, n_neighbors = 30
23:55:34 Building Annoy index with metric = cosine, n_trees = 50
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
23:55:40 Writing NN index file to temp file /tmp/RtmpyIf2G8/file2bdbd3ecb81da
23:55:40 Searching Annoy index using 1 thread, search_k = 3000
23:55:53 Annoy recall = 100%
23:55:55 Commencing smooth kNN distance calibration using 1 thread
23:55:59 Initializing from normalized Laplacian + noise
23:56:03 Commencing optimization for 200 epochs, with 1777116 positive edges
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
23:56:53 Optimization finished

[1] "allvar"
Start to do clustering for each sample...
>> Redoing truncated PCA on sample axolotl...
found nearest neighbors.
revoke Seurat to compute SNN.
done.
done. returning result...
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 16415
Number of edges: 645110

Running Louvain algorithm...

0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

Maximum modularity in 10 random starts: 0.9265
Number of communities: 21
Elapsed time: 2 seconds
>> Done clustering of sample axolotl.
>> Redoing truncated PCA on sample lizard...
found nearest neighbors.
revoke Seurat to compute SNN.
done.
done. returning result...
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 4187
Number of edges: 165059

Running Louvain algorithm...

0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

Maximum modularity in 10 random starts: 0.8960
Number of communities: 11
Elapsed time: 0 seconds
>> Done clustering of sample lizard.
>> Redoing truncated PCA on sample paintedturtle...
found nearest neighbors.
revoke Seurat to compute SNN.
done.
done. returning result...
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 18664
Number of edges: 622418

Running Louvain algorithm...

0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

Maximum modularity in 10 random starts: 0.8864
Number of communities: 12
Elapsed time: 2 seconds
>> Done clustering of sample paintedturtle.
Finished clustering.
Obtained average profiles of clusters.
Start to calculate standardized similarities to clusters...
Doing z-transformation...
Done.

00:04:34 UMAP embedding parameters a = 0.9922 b = 1.112
00:04:34 Read 39266 rows and found 44 numeric columns
00:04:34 Using Annoy for neighbor search, n_neighbors = 30
00:04:34 Building Annoy index with metric = cosine, n_trees = 50
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
00:04:40 Writing NN index file to temp file /tmp/RtmpyIf2G8/file2bdbd505d4a3f
00:04:40 Searching Annoy index using 1 thread, search_k = 3000
00:04:54 Annoy recall = 100%
00:04:55 Commencing smooth kNN distance calibration using 1 thread
00:04:59 Initializing from normalized Laplacian + noise
00:05:03 Commencing optimization for 200 epochs, with 1777102 positive edges
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
00:05:53 Optimization finished

[1] "commonvar"
Start to do clustering for each sample...
>> Redoing truncated PCA on sample axolotl...
found nearest neighbors.
revoke Seurat to compute SNN.
done.
done. returning result...
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 16415
Number of edges: 623330

Running Louvain algorithm...

0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

Maximum modularity in 10 random starts: 0.9174
Number of communities: 17
Elapsed time: 1 seconds
>> Done clustering of sample axolotl.
>> Redoing truncated PCA on sample lizard...
found nearest neighbors.
revoke Seurat to compute SNN.
done.
done. returning result...
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 4187
Number of edges: 172270

Running Louvain algorithm...

0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

Maximum modularity in 10 random starts: 0.8892
Number of communities: 11
Elapsed time: 0 seconds
>> Done clustering of sample lizard.
>> Redoing truncated PCA on sample paintedturtle...
found nearest neighbors.
revoke Seurat to compute SNN.
done.
done. returning result...
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 18664
Number of edges: 638406

Running Louvain algorithm...

0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

Maximum modularity in 10 random starts: 0.8901
Number of communities: 13
Elapsed time: 2 seconds
>> Done clustering of sample paintedturtle.
Finished clustering.
Obtained average profiles of clusters.
Start to calculate standardized similarities to clusters...
Doing z-transformation...
Done.

00:07:52 UMAP embedding parameters a = 0.9922 b = 1.112
00:07:52 Read 39266 rows and found 41 numeric columns
00:07:52 Using Annoy for neighbor search, n_neighbors = 30
00:07:52 Building Annoy index with metric = cosine, n_trees = 50
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
00:07:58 Writing NN index file to temp file /tmp/RtmpyIf2G8/file2bdbd60b3d33d
00:07:58 Searching Annoy index using 1 thread, search_k = 3000
00:08:13 Annoy recall = 100%
00:08:14 Commencing smooth kNN distance calibration using 1 thread
00:08:18 Initializing from normalized Laplacian + noise
00:08:21 Commencing optimization for 200 epochs, with 1810424 positive edges
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
00:09:11 Optimization finished

rPCA + Harmony, one2one ortholog

srat_rpca_harm_one_l = list()
for(n in names(srat_rpca_one_l)){
  print(n)
  # Run Harmony
  cell_integr = RunHarmony(srat_rpca_one_l[[n]], "species", tau = 30, 
                           plot_convergence = F, assay.use = "SCT")
  
  # Run UMAP on Harmony dimensions
  cell_integr = RunUMAP(cell_integr, reduction = "harmony", dims = 1:10)
  
  srat_rpca_harm_one_l[[n]] = cell_integr
}

CSS + Harmony, all ortholog

srat_css_harm_all_l = list()
for(n in names(srat_css_all_l)){
  print(n)
  # Run Harmony
  cell_integr = RunHarmony(srat_css_all_l[[n]], "species", tau = 30, reduction = "css",
                           dims.use = 1:ncol(Embeddings(srat_css_all_l[[n]],"css")),
                           plot_convergence = F, assay.use = "SCT")
  
  # Run UMAP on Harmony dimensions
  cell_integr = RunUMAP(cell_integr, reduction = "harmony", dims = 1:10)
  
  srat_css_harm_all_l[[n]] = cell_integr
}

CSS + Harmony, one2one ortholog

srat_css_harm_all_l = list()
for(n in names(srat_css_all_l)){
  print(n)
  # Run Harmony
  cell_integr = RunHarmony(srat_css_all_l[[n]], "species", tau = 30, reduction = "css",
                           dims.use = 1:ncol(Embeddings(srat_css_all_l[[n]],"css")),
                           plot_convergence = F, assay.use = "SCT")
  
  # Run UMAP on Harmony dimensions
  cell_integr = RunUMAP(cell_integr, reduction = "harmony", dims = 1:10)
  
  srat_css_harm_all_l[[n]] = cell_integr
}

[1] "all"

Harmony 1/10
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 2/10
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 3/10
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 4/10
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 5/10
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 6/10
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 7/10
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony converged after 7 iterations
Invalid name supplied, making object name syntactically valid. New object name is Seurat..ProjectDim.SCT.harmony; see ?make.names for more details on syntax validity00:11:22 UMAP embedding parameters a = 0.9922 b = 1.112
00:11:22 Read 39266 rows and found 10 numeric columns
00:11:22 Using Annoy for neighbor search, n_neighbors = 30
00:11:22 Building Annoy index with metric = cosine, n_trees = 50
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
00:11:27 Writing NN index file to temp file /tmp/RtmpyIf2G8/file2bdbd3ea2902b
00:11:27 Searching Annoy index using 1 thread, search_k = 3000
00:11:43 Annoy recall = 100%
00:11:45 Commencing smooth kNN distance calibration using 1 thread
00:11:48 Initializing from normalized Laplacian + noise
00:11:50 Commencing optimization for 200 epochs, with 1651614 positive edges
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
00:12:39 Optimization finished

[1] "allvar"

Quick-TRANSfer stage steps exceeded maximum (= 1963300)Harmony 1/10
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 2/10
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 3/10
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 4/10
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 5/10
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 6/10
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 7/10
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 8/10
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 9/10
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 10/10
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Invalid name supplied, making object name syntactically valid. New object name is Seurat..ProjectDim.SCT.harmony; see ?make.names for more details on syntax validity00:14:46 UMAP embedding parameters a = 0.9922 b = 1.112
00:14:46 Read 39266 rows and found 10 numeric columns
00:14:46 Using Annoy for neighbor search, n_neighbors = 30
00:14:46 Building Annoy index with metric = cosine, n_trees = 50
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
00:14:51 Writing NN index file to temp file /tmp/RtmpyIf2G8/file2bdbd41e214cd
00:14:51 Searching Annoy index using 1 thread, search_k = 3000
00:15:06 Annoy recall = 100%
00:15:07 Commencing smooth kNN distance calibration using 1 thread
00:15:10 Initializing from normalized Laplacian + noise
00:15:12 Commencing optimization for 200 epochs, with 1645504 positive edges
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
00:16:01 Optimization finished

[1] "commonvar"

Harmony 1/10
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 2/10
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 3/10
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 4/10
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 5/10
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 6/10
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 7/10
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony converged after 7 iterations
Invalid name supplied, making object name syntactically valid. New object name is Seurat..ProjectDim.SCT.harmony; see ?make.names for more details on syntax validity00:17:37 UMAP embedding parameters a = 0.9922 b = 1.112
00:17:37 Read 39266 rows and found 10 numeric columns
00:17:37 Using Annoy for neighbor search, n_neighbors = 30
00:17:37 Building Annoy index with metric = cosine, n_trees = 50
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
00:17:42 Writing NN index file to temp file /tmp/RtmpyIf2G8/file2bdbd5e63e6fa
00:17:42 Searching Annoy index using 1 thread, search_k = 3000
00:17:58 Annoy recall = 100%
00:17:59 Commencing smooth kNN distance calibration using 1 thread
00:18:03 Initializing from normalized Laplacian + noise
00:18:05 Commencing optimization for 200 epochs, with 1670228 positive edges
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
00:18:54 Optimization finished

Save lists

Save

srat_css_harm_one_l = list()
for(n in names(srat_css_one_l)){
  print(n)
  # Run Harmony
  cell_integr = RunHarmony(srat_css_one_l[[n]], "species", tau = 30, reduction = "css",
                           dims.use = 1:ncol(Embeddings(srat_css_one_l[[n]],"css")),
                           plot_convergence = F, assay.use = "SCT")
  
  # Run UMAP on Harmony dimensions
  cell_integr = RunUMAP(cell_integr, reduction = "harmony", dims = 1:10)
  
  srat_css_harm_one_l[[n]] = cell_integr
}

[1] "all"

Harmony 1/10
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 2/10
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 3/10
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 4/10
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 5/10
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony converged after 5 iterations
Invalid name supplied, making object name syntactically valid. New object name is Seurat..ProjectDim.SCT.harmony; see ?make.names for more details on syntax validity00:20:34 UMAP embedding parameters a = 0.9922 b = 1.112
00:20:34 Read 39266 rows and found 10 numeric columns
00:20:34 Using Annoy for neighbor search, n_neighbors = 30
00:20:34 Building Annoy index with metric = cosine, n_trees = 50
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
00:20:40 Writing NN index file to temp file /tmp/RtmpyIf2G8/file2bdbd293b9c4c
00:20:40 Searching Annoy index using 1 thread, search_k = 3000
00:20:53 Annoy recall = 100%
00:20:55 Commencing smooth kNN distance calibration using 1 thread
00:20:58 Initializing from normalized Laplacian + noise
00:21:00 Commencing optimization for 200 epochs, with 1654888 positive edges
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
00:21:49 Optimization finished

[1] "allvar"

Harmony 1/10
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 2/10
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 3/10
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 4/10
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 5/10
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 6/10
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 7/10
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 8/10
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 9/10
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony converged after 9 iterations
Invalid name supplied, making object name syntactically valid. New object name is Seurat..ProjectDim.SCT.harmony; see ?make.names for more details on syntax validity00:23:51 UMAP embedding parameters a = 0.9922 b = 1.112
00:23:51 Read 39266 rows and found 10 numeric columns
00:23:51 Using Annoy for neighbor search, n_neighbors = 30
00:23:51 Building Annoy index with metric = cosine, n_trees = 50
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
00:23:56 Writing NN index file to temp file /tmp/RtmpyIf2G8/file2bdbd4c0ca239
00:23:56 Searching Annoy index using 1 thread, search_k = 3000
00:24:10 Annoy recall = 100%
00:24:11 Commencing smooth kNN distance calibration using 1 thread
00:24:14 Initializing from normalized Laplacian + noise
00:24:16 Commencing optimization for 200 epochs, with 1651456 positive edges
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
00:25:05 Optimization finished

[1] "commonvar"

did not converge in 25 iterationsHarmony 1/10
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 2/10
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 3/10
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 4/10
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 5/10
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony converged after 5 iterations
Invalid name supplied, making object name syntactically valid. New object name is Seurat..ProjectDim.SCT.harmony; see ?make.names for more details on syntax validity00:26:31 UMAP embedding parameters a = 0.9922 b = 1.112
00:26:31 Read 39266 rows and found 10 numeric columns
00:26:31 Using Annoy for neighbor search, n_neighbors = 30
00:26:31 Building Annoy index with metric = cosine, n_trees = 50
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
00:26:36 Writing NN index file to temp file /tmp/RtmpyIf2G8/file2bdbd88888be
00:26:36 Searching Annoy index using 1 thread, search_k = 3000
00:26:51 Annoy recall = 100%
00:26:52 Commencing smooth kNN distance calibration using 1 thread
00:26:55 Initializing from normalized Laplacian + noise
00:26:57 Commencing optimization for 200 epochs, with 1660784 positive edges
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
00:27:46 Optimization finished

Plots for all genes

#saveRDS(srat_rpca_all_l, file = "./data/integration/srat_rpca_all_l.RDS")
#saveRDS(srat_rpca_one_l, file = "./data/integration/srat_rpca_one_l.RDS")

#saveRDS(srat_harm_all_l, file = "./data/integration/srat_harm_all_l.RDS")
#saveRDS(srat_harm_one_l, file = "./data/integration/srat_harm_one_l.RDS")

saveRDS(srat_css_all_l, file = "./data/integration/srat_css_all_l.RDS")
saveRDS(srat_css_one_l, file = "./data/integration/srat_css_one_l.RDS")

#saveRDS(srat_rpca_harm_all_l, file = "./data/integration/srat_rpca_harm_all_l.RDS")
#saveRDS(srat_rpca_harm_one_l, file = "./data/integration/srat_rpca_harm_one_l.RDS")

saveRDS(srat_css_harm_all_l, file = "./data/integration/srat_css_harm_all_l.RDS")
saveRDS(srat_css_harm_one_l, file = "./data/integration/srat_css_harm_one_l.RDS")

Plots for one2one genes

for(n in names(srat_rpca_all_l)){
  plt1 = DimPlot(srat_rpca_all_l[[n]], reduction = "umap", group.by = "species")+
    labs(title = paste0("rPCA, all orthologs, ", n))+
  DimPlot(srat_rpca_all_l[[n]], reduction = "umap", group.by = "classes")+
    labs(title = paste0("rPCA, all orthologs, ", n))
  print(plt1)
  
  plt2 = DimPlot(srat_css_all_l[[n]], reduction = "umap_css", group.by = "species")+
    labs(title = paste0("CSS, all orthologs, ", n))+
  DimPlot(srat_css_all_l[[n]], reduction = "umap_css", group.by = "classes")+
    labs(title = paste0("CSS, all orthologs, ", n))
  print(plt2)
  
  plt3 = DimPlot(srat_harm_all_l[[n]], reduction = "umap", group.by = "species")+
    labs(title = paste0("Harmony, all orthologs, ", n))+
  DimPlot(srat_harm_all_l[[n]], reduction = "umap", group.by = "classes")+
    labs(title = paste0("Harmony, all orthologs, ", n))
  print(plt3)
  
  plt4 = DimPlot(srat_rpca_harm_all_l[[n]], reduction = "umap", group.by = "species")+
    labs(title = paste0("rPCA+Harmony, all orthologs, ", n))+
  DimPlot(srat_rpca_harm_all_l[[n]], reduction = "umap", group.by = "classes")+
    labs(title = paste0("rPCA+Harmony, all orthologs, ", n))
  print(plt4)
  
  plt5 = DimPlot(srat_css_harm_all_l[[n]], reduction = "umap", group.by = "species")+
    labs(title = paste0("CSS+Harmony, all orthologs, ", n))+
  DimPlot(srat_css_harm_all_l[[n]], reduction = "umap", group.by = "classes")+
    labs(title = paste0("CSS+Harmony, all orthologs, ", n))
  print(plt5)
}